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ABSTRACT 

We present distributions of the orbital parameters of dark matter substructures at 
the time of merging into their host halo. Accurate knowledge of the orbits of dark 
matter substructures is a crucial input to studies which aim to assess the effects of the 
cluster environment on galaxies, the heating of galaxy disks and many other topics. 
Orbits are measured for satellites in a large number of N-body simulations. We focus 
on the distribution of radial and tangential velocities, but consider also distributions 
of orbital eccentricity and semi-major axis. We show that the distribution of radial 
and tangential velocities has a simple form and provide a fitting formula for this 
distribution. We also search for possible correlations between the infall directions of 
pairs of satellites, finding evidence for positive correlation at small angular separations 
as expected if some infall occurs along filaments. We also find (weak) evidence for 
correlations between the direction of the infall and infall velocity and the spin of the 
host halo. 
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1 INTRODUCTION 



In currently favoured cosmological models, dark matter ha- 
los grow via the merging together of smaller systems, lead- 
ing to an ever-growing hierarchy of halos. Recent numerical 
simulations have demonstrated that the remnants of pre- 
existing dark matter halos which merged to become part of 
a larger system (the "host") can survive for significant pe- 
riods of time within the larger system (Moore et al. 1999; 
Klypin et al. 1998) . These subhalos orbit around in the po- 
tential of the host gradually losing mass via tidal forces and 
spiralling in to ever smaller radii due to dynamical friction. 
These substructures (or at least some subset of them) are 
presumably the abodes of satellite galaxies, such as those 
found in the Local Group, and of the majority of cluster 
galaxies. 

This substructure has attracted a great deal of inter- 
est since its discovery. Observational tests for its presence, 
though not yet conclusive, are in good agreement with the 
theoretical expectations (Metcalf & Madau 2001; Chiba 
2002; Dalai & Kochanek 2002). There has been much work 
conducted in which the distribution and properties of sub- 
structures, their effects on galaxy disks and so on were exam- 
ined (Ghigna et al. 1998; Tormen, Diaferio & Syer 1998; van 
den Bosch et al. 1999; Zhang et al. 2002; Benson et al. 2004; 
Gao et al. 2004; Diemand, Moore & Stadel 2004). While 
the orbital parameters of substructures have been measured 
in the past this measurement has often been at the end point 
of the substructure evolution (i.e. at the present day) when 
significant dynamical evolution in the orbital parameters is 



expected (e.g. Ghigna et al. 1998). Exceptions to this are 
the works of Tormen (1997), Vitvitska et al. (2002) and 
Khochfar & Burkert (2004). Tormen (1997) and Khoch- 
far & Burkert (2004) both identified progenitors of halos 
in their N-body simulations and measured the orbital pa- 
rameters of them, while Vitvitska et al. (2002) searched for 
pairs of halos about to merge and measured the orbital pa- 
rameters of these. These works have typically made use of 
rather small samples of orbits and (perhaps consequently) 
have been unable to fully characterise the two dimensional 
distribution of orbital parameters needed to construct real- 
istic orbits. 



The distribution of initial orbital parameters of sub- 
structure halos at the time of merging into the host sys- 
tem is a particularly interesting property as it represents 
the initial conditions which determine the later evolution of 
the substructure within the host. The effectiveness of many 
processes invoked to explain the morphological transforma- 
tion of galaxies in clusters (e.g. ram pressure stripping, tidal 
mass loss, galaxy harassment, etc.) depend crucially on the 
nature of the galaxy orbit (see, for example, Moore, Lake & 
Katz 1998; Abadi, Bower & Navarro 2000). The distribu- 
tion of orbits will also determine the rate of galaxy mergers 
and therefore the degree of heating and rate of morpholog- 
ical transformation experienced by galaxy disks. Taking a 
more practical point of view, recent semi-analytic models of 
satellite halo orbits (Benson et al. 2002; Taylor & Babul 
2004) have been able to follow the orbital evolution of satel- 
lites quite accurately, but these models are only as good as 
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their initial conditions which, until now, have been known 
only rather poorly. 

In this work, we quantify the distribution of orbital pa- 
rameters for dark matter halos at the point of merging with 
their host (i.e. we proceed in a similar way as did Vitvitska et 
al. 2002). We measure this distribution in a large number of 
N-body simulations to attain high statistical precision and 
to facilitate checks of our techniques and tests for variations 
of the distribution of orbital parameters with variables such 
as redshift, halo mass etc. While we will present distribu- 
tions of orbital eccentricity and semi-major axis, our focus 
is on distributions of radial and tangential velocities, which 
we find are more practical when dealing with orbits in non- 
spherical systems in which dynamical friction is at work 1 . 
We also examine the distribution of infalling substructures 
as a function of position on the virial sphere, and explore 
correlations between orbital properties and the spin of the 
host halo. 

Our aim is to provide a precise and accurate measure- 
ment of the distribution of orbital properties of substruc- 
tures at the time of merging, and to provide fits to this 
distribution so that it may be used in further studies. This 
distribution could, in principle, depend on many quantities, 
such as the masses of the merging halos, redshift, cosmologi- 
cal parameters etc. Furthermore, the six parameters describ- 
ing each orbit (e.g. the position and velocity of the satellite 
at the time of merging, or any equivalent parameter set) 
may well be correlated with each other, such that we should 
really examine a six-dimensional phase-space distribution 
function. With the currently available N-body simulations 
we will limit ourselves to exploring a two-dimensional func- 
tion, typically that of radial and tangential velocities (ef- 
fectively assuming that infalling satellites are uniformly dis- 
tributed on a sphere around the halo centre and that their 
tangential velocities have no preferred direction), although 
we will explore correlations between these quantities and the 
host halo. We note also that the situation could in principle 
be more complicated still. We are aiming to quantify P(x), 
where x are the orbital parameters and P is the distribu- 
tion of these averaged over all merging events. However, af- 
ter one merger with parameters xi the relevant distribution 
function for the next merger may be different, P(xjxi). An 
example might be infall of halos along a filament. Knowing 
that one halo fell in from a particular direction, it becomes 
more likely that the next halo will fall in from a similar 
direction. We will explore one aspect of this possibility by 
measuring the distribution of angles between pairs of in- 
falling satellites. 

The remainder of this paper is arranged as follows. In 
33 we describe our analysis technique while in JHlwe present 
our results. We give our conclusions in 



2 ANALYSIS 

2.1 N-body Simulations 

To measure satellite orbital parameters we make use of 
a large number of N-body simulations carried out by the 
VIRGO Consortium and which are publicly available (see 
Jenkins et al. 1998; Kauffmann et al. 1999.; Jenkins et al. 
2001 for further details), together with one other simulation 
used for testing various aspects of our methodology. These 
span a range of cosmologies and redshifts. Details of the sim- 
ulations used are given in Table Q All of the outputs from 
these simulations are analysed, but in practise only those at 
redshifts z <J2 provide statistically useful measurements of 
orbital parameter distributions. 



2.2 Group Finding 

In order to find merging dark matter halos in the simula- 
tions we must first identify all dark matter halos. To locate 
dark matter halos in the N-body simulations we employ two 
standard group finders, the friends-of-friends (FOF; Davis 
et al. 1985) and spherical overdensity (SO; Lacey & Cole 
1994) algorithms. We will compare results for halos found 
using these two techniques to test for any dependence on 
the group finding algorithm used. 

Each algorithm has one tunable parameter, the linking 
length, rii n k, for the FOF algorithm and the mean density 
contrast inside the sphere, A, for the SO algorithm. Both can 
be related to the mean density of dark matter halos (once 
a specific halo density profile has been chosen in the case of 
the FOF algorithm). We apply each algorithm twice, once 
assuming a mean overdensity for halos of A = 187r 2 ~ 177.7 
(equivalent to riink — 0.20f, assuming an isothermal halo 
profile 2 , where f is the mean inter-particle spacing in the 
simulation), as expected from the spherical collapse model 
in a critical density cosmology (e.g. Peebles 1980), and once 
using the mean overdensity expected from the spherical col- 
lapse model for the specific cosmology and redshift in ques- 
tion (Lacey & Cole 1993; Eke, Cole & Frenk 1996). We will 
refer to these two alternatives as "fixed A" and "variable A" 
respectively, and will compare results from the two. 

Once halos have been located by either algorithm we 
apply a procedure to remove unbound halos from the re- 
sulting catalogue. Our technique is described fully by Ben- 
son et al. (2001) and involves repeatedly removing the least 
bound particle from an unbound halo until the halo either 
becomes bound, or falls below the minimum mass required 
to be included in our catalogue. 



1 Since the orbital parameters are constantly changing for such 
orbits, the eccentricity and peri-centric distance no longer have 
the advantage of being constant along the orbit. The orbital ve- 
locities are more closely related to the quantities required by 
semi-analytic orbital models so we prefer to use them. The two 
pairs of orbital parameters (eccentricity+semi-major axis and ra- 
dial+tangential velocity) are, of course, equivalent. 



2 It is well known that cold dark matter halos are not well ap- 
proximated by isothermal spheres. However, if we instead adopt 
an NFW density profile (Navarro, Frenk & White 1997) for our 
halos the appropriate value of ri; n k ranges between 0.22F and 
0.26f for halos with concentrations in the range 5 to 15. As such, 
a somewhat larger value of may be appropriate. Neverthe- 

less, we will retain the convention of assuming isothermal halos 
here and resign a study of the most appropriate linking length to 
use to future work. 
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Table 1. The names, parameters and output redshifts of the N-body simulations used in our analysis. The first two columns give the 
name of the simulation set and the cosmological model respectively. Columns 3 lists the number of particles in each simulation, while 
columns 4 and 5 list the cosmological parameters Qo and Ao appropriate to each simulation. Column 6 specifies the length of the 
simulation cube, while column 7 specifies the mass of each particle in the simulation. Column 8 gives the softening length used in the 
simulation. Finally, column 9 lists the redshifts at which outputs from the simulation are available. 



Simulation Model Particles Qq Aq L/h x Mpc m p /h 1 Mq l ao {th x kpc Redshifts 



GIF 


ACDM 


256 3 


0.3 


0.7 


141.3 


1.4 x 10 10 


20 


50, uniform in ln(a) from z = 50 to z = 


GIF 


OCDM 


256 3 


0.3 


0.0 


141.3 


1.4 x 10 10 


30 


0.0, 0.1, 0.3, 0.5, 1.0, 1.5, 2.0, 3.0, 5.0 


GIF 


SCDM 


256 3 


1.0 


0.0 


84.5 


1.0 x 10 10 


36 


0.0, 0.1, 0.3, 0.5, 1.0, 1.5, 2.0, 3.0, 5.0 


GIF 


tCDM 


256 3 


1.0 


0.0 


84.5 


1.0 x 10 10 


36 


0.0, 0.1, 0.3, 0.5, 1.0, 1.5, 2.0, 3.0, 5.0 


GIF-n 


rCDM 


256 3 


1.0 


0.0 


84.5 


1.0 x 10 10 


36 


0.0 


Virgo 


ACDM 


256 3 


0.3 


0.7 


239.5 


6.86 x 10 10 


25 


0.0, 0.1, 0.3, 0.5, 1.0, 1.5, 2.0, 3.0, 5.0 


Virgo 


OCDM 


256 3 


0.3 


0.0 


239.5 


6.85 x 10 10 


30 


0.0, 0.1, 0.3, 0.5, 1.0, 1.5, 2.0, 3.0, 5.0 


Virgo 


SCDM 


256 3 


1.0 


0.0 


239.5 


2.27 x 10 11 


36 


0.0, 0.1, 0.3, 0.5, 1.0, 1.5, 2.0, 3.0, 5.0 


Virgo 


rCDM 


256 3 


1.0 


0.0 


239.5 


2.27 x 10 11 


36 


0.0, 0.1, 0.3, 0.5, 1.0, 1.5, 2.0, 3.0, 5.0 


ff VLS 


ACDM 


512 3 


0.3 


0.7 


479.0 


6.86 x 10 10 


30 


0.0, 0.5, 1.0, 2.0, 3.0, 5.0 



2.3 Defining the Halo Centre 

To measure orbital properties of infalling satellites we need 
to define a centre (both in position and velocity) for each 
halo in order to have a suitable origin for our coordinate 
system. The simplest option is to determine the centre of 
mass and the mass weighted mean velocity of the halo and 
take these as the origin. We call this "COM centring". Be- 
cause of its simplicity we will examine results based upon 
this approach. However, while a simple centre of mass es- 
timate of the halo centre is reasonable if halos are smooth, 
spherical systems, in reality it has many failings (particu- 
larly when applied to FOF halos) . The FOF algorithm often 
links together halos that are about to merge by a low den- 
sity "bridge" of particles. This will skew the centre of mass 
of the halo away from what perhaps should be considered 
the centre (e.g. the position corresponding to the centre of 
mass of the main component of the merging system). Be- 
cause of this limitation we will adopt a second approach in 
which we define the centre of a halo as being the position 
of the particle with the lowest gravitational energy (count- 
ing only interactions with other particles in the halo). This 
will naturally pick out a particle in the densest region and, 
given two halos joined by a low-density bridge should pick 
out a particle in the more massive of the two halos. However, 
we cannot take the velocity of this particle as being repre- 
sentative of the velocity of the halo, since its motion will 
consist of the mean halo motion plus a component due to 
the halo's internal velocity dispersion. Unfortunately, just 
as in position space, halos in velocity space can show bi- 
modal distributions (as happens when a halo is linked to 
a nearby halo which is infalling). This can bias the mass 
weighted mean velocity estimate of the origin away from 
the "correct" value. To circumvent this problem we adopt 
a similar approach in velocity space as in position space. 
Namely, we estimate an "energy", e, for each particle, i, us- 
ing ti — V* . ,. — 1/lvj — Vjl, with the sum taken over all 
particles in the halo, and then locate the particle with the 
lowest energy. This should lie close to the true mean motion 
of the halo. We call this method "MBP centring". 

It is worth noting that the velocity origin can differ 
significantly between the two definitions we adopt. Fig. Q 
shows the centre of a particular halo from the z = output 
of the GIF ACDM simulation in both position and velocity 



space. Each frame has its origin on the most-bound parti- 
cle, as marked by the dashed lines, while the dotted lines 
indicate the centre of mass or mass- weighted mean velocity. 
In this example, the centre of mass almost coincides with 
the most bound particle (somewhat fortuitously as nearby 
halos on either side, linked in by the FOF algorithm, are 
cancelling each other out). In velocity space however, we see 
that the velocity origin is shifted by over 500km/s from the 
more realistic velocity origin. This could seriously affect our 
estimates of orbital parameters. 



2.4 Satellite Orbital Parameters 

From our catalogue of dark matter halos in each simulation 
we search for pairs of halos which are about to merge. From 
here on, all velocities are measured in units of the virial 
circular velocity of the host halo, V v ir, and all radii in units 
of virial radius of the host halo, r v i r , as we expect these to 
be characteristic scales of the systems 3 . We search for halos 
within a distance from the host centre between r = 1 ± Ar, 
and which have an inward directed velocity, v (i.e. r ■ v < 
where r is the vector from the centre of the host to the centre 
of the potential satellite halo). We choose Ar = 0.2. Note 
that we allow for the possibility of halos with r < 1 since 
the non-spherical shape of real halos can permit a halo to 
remain separate from the host even when r < 1. It should 
be noted that this radial selection biases us against finding 
mergers between halos of comparable mass (since in this case 
it is unlikely that the satellite will remain as an isolated halo 
once its centre is within 1 + Ar). For present purposes this 
bias is unimportant, and so we retain the above criterion for 
simplicity. This bias could however, be easily circumvented 
by adopting a radial selection based upon the sum of the host 
and satellite virial radii instead. From the halos selected in 
this way, we compute the radial and tangential components 



3 We convert the comoving coordinates of the N-body simulation 
to physical coordinates and add on the Hubble flow to the peculiar 
velocities taken from the N-body simulations. Halo virial radii and 
velocities are determined from their masses assuming the halo to 
have the mean density appropriate to a just-collapsed spherical 
top-hat overdensity. 
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Figure 1. The results of using the COM and MBP algorithms to define the origin of the coordinate system in a dark matter halo 
identified in the z = output of the GIF ACDM simulation. The upper row shows three projections of the spatial distribution of 
particles. The intersection of the dashed lines indicates the origin according to the MBP algorithm, while that of the dotted line indicates 
the origin according to the COM algorithm. The lower row shows projections of the same halo in velocity space. Dashed and dotted lines 
are as in the upper row. 



of velocity. We also store the three dimensional position and 
angular momentum of the merging satellite. 

Since we are interested in the orbital parameters of 
satellites as they cross the virial radius of a larger halo we 
correct our orbital parameters (which are measured at some 
radius close to, but not equal to, the virial radius). To do this 
we treat the two halos as point masses, and simply determine 
the point at which the satellite's orbit first crosses the virial 
radius of the larger halo. We store the position, velocity and 
angular momentum of the satellite at this point. This ap- 
proach is an approximation for two reasons. Firstly, as the 
host halo is not a point mass, the mass interior to the sub- 
structure's orbit will change along that orbit. In practise this 
effect is quite small, leading to only a 5% error in the orbital 
velocities. (Note also that the density profile is not spheri- 
cally symmetric, which will lead to further errors.) Secondly, 
we neglect the effects of dynamical friction on the orbital 
parameters. A simple estimate based upon Chandrasekhar's 
methods indicates that this would lead to an error in our 
orbital velocities of around 10% for a substructure to host 



mass ratio of 0.08 (which is typical of the systems found 
in our simulations), scaling approximately in proportion to 
this ratio. All of these problems could be largely overcome 
by solving the equations of motion for the substructure in a 
realistic host potential including a dynamical friction term. 
This will be the focus of future work. 

Some fraction of substructures are found to be on un- 
bound orbits. This presents no problem for our analysisb, 
we can of course still measure the orbital parameters of such 
substructures, and so we retain these objects in our calcula- 
tions. The fate of such substructures will be discussed below. 
Some substructures are found with r < 1 — already inside 
the virial radius by our definition. These substructures are 
propagated backwards along their orbit to find their orbital 
parameters at the time of crossing r — 1 (as with all or- 
bits, no account is made for any mass loss which might have 
occurred from these halos, nor for the effects of dynamical 
friction). Finally, we find some halos whose orbits do not 
cross the virial radius of the host. Such halos are flagged 
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as being "bad" and are treated separately from other halos 
(see jnj. 

We must also account for the fact that our selection of 
halos with 1 — Ar < r < 1 + Ar leads to a bias against finding 
radial orbits as they will spend less time in this region than 
more circular orbits. To correct for this we simply deter- 
mine, from the measured orbital parameters of the satellite, 
the time, St, it takes to traverse the region r — 1+Ar to r m i u . 
Here r m ] n is the minimum radius at which the satellite halo 
would have been identified by the group finder. When con- 
structing distributions of orbital parameters we then weight 
by At /St where At is the cosmic time between the current 
N-body simulation output and the previous one (or t — in 
the case of the highest redshift output). 

The determination of r m i n depends on the group finder 
used. With the SO group finder it is relatively easy to de- 
termine r m i n . Under the SO algorithm each halo is assigned 
a radius (the radius containing a mean overdensity of some 
specified value). Once all halos have been found any halo 
whose centre lies within the radius, rso, of a larger halo is 
merged with that larger halo and removed from the list of 
individual halos. (Note that the radius of the larger halo is 
not changed by this merging.) Thus, r m j n is simply rso, or 
1 — Ar, whichever is larger. 

For the FOF group finder things are a little more com- 
plicated as the halos found are not spherical. The satellite 
halo would no longer have been found as an isolated object 
by the group finding algorithm once any one of its parti- 
cles came within a distance rii n k of a particle in the larger 
halo. We therefore search for the first point along the or- 
bit of the satellite at which any one of its particles comes 
within rii n k of a particle in the larger halo. We define r m i n 
to be the orbital radius at this point, or 1 — Ar, whichever is 
larger. The advantage of this approach is that it works even 
for the non-spherical halos found by the FOF algorithm. Its 
disadvantage is that it treats the orbit as that of two point 
masses and also ignores any internal evolution of the satellite 
or host halos during the time it takes the satellite to move 
along its orbit. This latter is not a problem providing the 
two halos are in internal equilibrium and not rotating since 
then, although the individual particles in the halos move, 
their distribution at any time provides a fair sample of the 
mass distribution of the halo at any later time. Of course, 
in reality the halos will not be in equilibrium (although we 
expect them to be close to it). In particular, the FOF algo- 
rithm is known to make "dumbbell-shaped" halos by linking 
together two halos by a low density bridge. These are cer- 
tainly not equilibrium systems in the sense used here. They 
are also those in which the two-body orbit approximation is 
likely to be worst. We consider this to be a limitation of the 
FOF algorithm, and do not explore more complicated ways 
of dealing with this problem here. 

It should be noted that, with our method for locating 
merging halos, some host halos may be experiencing mergers 
with multiple substructrues at any given time. In fact, we 
find that about 25% of all of our merger events at z = in- 
volve two or more substructures accreting onto the same 
host halo. For the largest clusters we find up to around 
twenty ongoing mergers in some cases. We find very few 
mergers with low mass ratios (e.g. less than 4:1). As such, 
the inclusion or not of hosts currently underdoing major 
mergers does not affect our results significantly. 



3 RESULTS 

We examine the orbital parameter distributions for each in- 
dividual output of each simulation. We will also combine 
results together where possible to improve the statistical 
precision. All results will make use of the FOF halo finding 
algorithm, MBP halo centring and the variable A method 
for setting rrmk/A unless otherwise stated. 

Figure |5] shows an example of the distribution of or- 
bital parameters that we measure. The distribution of radial 
and tangential velocities (upper left and right-hand panels 
respectively) have quite simple, and perhaps unsurprising, 
forms, being peaked at V ~ 1 with a dispersion of order 
unity. The infall angle, defined as the (negative of the) an- 
gle between the infalling substructure's radius and velocity 
vectors (i.e. <f> = — cos -1 r • v/|r|/|v|), is shown in the lower- 
left hand panel. This distribution will be investigated fur- 
ther in H3.3I Finally, the lower right-hand panel shows the 
two-dimensional distribution of radial and tangential orbital 
velocities. It is clear that there is a significant correlation be- 
tween these two parameters. Another interesting feature of 
this distribution is that a significant fraction of orbits drawn 
from this distribution are initially unbound. The energy of 
orbits, in our units, is given by 

E - l+ wX v * +{ ^ff vi )> (1) 

where f 2 = 1 + M 2 /M 1 . Note that f 2 = M 2 /fi where fj, = 
M 1 M 2 /(M 1 + M 2 ) is the usual reduced mass. The dotted 
line in Fig. |5] shows the line E — for the case f 2 = l (i.e. 
Mi <C M 2 ). Points to the upper right of this line correspond 
to unbound orbits. For the particular distribution shown 
about 18% of all orbits are unbound. We choose to retain 
these orbits for two reasons: 

(i) When using the measured distribution to select initial 
orbits for satellites, unbound orbits can easily be discarded 
if desired. 

(ii) Due to the effects of dynamical friction, an orbit that 
starts out unbound will not necessarily stay that way. 

To examine the importance of this second point we employ 
the semi-analytic model of Benson et al. (2004) which fol- 
lows the cosmological growth of dark matter halos (and their 
associated galaxies) including a detailed treatment of the 
orbital evolution of satellite halos. In Benson et al. (2004) 
the initial orbits of merging satellites were determined by 
setting the energy of each orbit equal to that of a circu- 
lar orbit at half the virial radius and choosing a circularity 
(i.e. the angular momentum of the satellite in units of that 
of a circular orbit with the same energy) from a uniform 
distribution between 0.1 and 1.0. These choices were mo- 
tivated by the results of Ghigna et al. (1998). Here, we 
instead use the measured distribution of orbital velocities, 
including unbound orbits, to set the initial velocity of satel- 
lites, and choose their initial position at random on a sphere 
with radius equal to the virial radius of their host. From 
this cosmologically representative sample of halos and or- 
bits, we identify those which start out unbound. Of these, 
some fraction will lose sufficient energy through dynamical 
friction that they become bound by the endpoint of their 
evolution (i.e. by z = 0) while others will fail to do so and 
will instead leave their host halo with positive energy. We 
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Figure 2. Distributions of orbital parameters measured in the 
VLS plus VIRGO ACDM z = output. Upper left and right- 
hand panels show distributions of radial and tangential veloci- 
ties respectively. The lower left-hand panel shows the distribu- 
tion of infall angles, while the lower right-hand panel shows the 
two-dimensional distribution of radial and tangential velocities. 
Contours are drawn at d 2 //dVrdV 9 = 0.01, 0.1, 0.5, 1.0 and 1.4 
from lightest to heaviest lines. The division between bound and 
unbound orbits in this panel is shown by the dotted line. 

find that approximately 2% of all initially unbound orbits 
(equivalent to 0.3% of all orbits) fail to become bound and 
so escape their halo. As such, these "lost" satellites are only 
a small fraction of the total. 

Our results are in good agreement with previous work. 
Figure |3] shows a comparison of the distribution of tangen- 
tial velocities with that found by Vitvitska et al. (2002; our 
Ve is equivalent to their L/L v i T ). Although the two distribu- 
tions differ as judged by a % 2 test, the discrepancy is due to 
two points and plausibly reflects our ignorance of the true 
errors and the differences in the simulations (e.g. softening, 
method of force calculation etc.) used in this work and that 
of Vitvitska et al. (2002). 

3.1 Tests of the Distributions 

Firstly, we examine which, if any of our measured distribu- 
tions are consistent with each other. This will allow us to 
determine which distributions we can realistically average 
together in order to improve the statistical precision of our 
measurements. 

3.1.1 Calibration of x 2 

We adopt a simple x 2 test to determine if two of our 
measured two-dimensional velocity distributions are consis- 
tent with each other. It should be noted that the errors 
which we determine for our distributions are likely to be an 
underestimate — they account for the finite number of merg- 
ers in each bin, but ignore such contributions as errors in 



Figure 3. The distribution of tangential velocities for orbits. Cir- 
cles show results for the VLS plus VIRGO ACDM simulations 
z = output from this work, while crosses (offset horizontally 
slightly for clarity) show the results of Vitvitska et al. (2002). 



our determinations of orbital velocities etc. Given this, and 
the fact that our errors may not be normally distributed, we 
would ideally like a calibration of the x 2 test. To achieve this 
we compare distributions from our GIF and GIF-n rCDM 
z = simulations. Comparing both the FOF and SO results, 
with halo centres defined using both centre of mass and most 
bound algorithms we find values of x 2 P er degree of freedom 
which scatter around unity, with a mean of 1.05. Although 
we would ideally like many more independent simulations to 
test our errors this gives us confidence that the errors are 
a good approximation to the true uncertainty on each data 
point. 



3.1.2 Distribution With and Without "Bad" Orbits 

A small fraction of the orbits that we find are flagged as 
being "bad" in the sense that they do not pass through one 
or both of the radial limits which we use for computing the 
weight to assign to each orbit. This may represent cases in 
which a halo formed within the outer radial limit (and so 
never passed through it), or, more likely, a limitation of the 
simple, two-body orbit neglecting dynamical friction that 
we use to approximate the motion of the halos. The best 
guess at a suitable weight for these orbits is to use their in- 
stantaneous radial velocity to determine the time taken to 
cross between the two radial limits. However, we find that 
the resulting distributions of orbital parameters for bad or- 
bits differ significantly (as judged by the \ 2 test) from those 
of good orbits. Therefore, we adopt the approach of excis- 
ing all bad orbits from our distributions. Ideally, we should 
deal with these better by solving for the orbit correctly (i.e. 
including extended masses and dynamical friction) to see if 
they really do merge and thereby assigning a realistic weight. 
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3.1.3 Number of particles per halo 

Our halo finding algorithms retain only halos consisting of 
ten particles or more. To test whether particle number has 
any effect on the measured distribution of orbit parame- 
ters we compare measurements of the orbit distribution in 
the VLS and VIRGO simulations with the equivalent GIF 
simulations, keeping halos with 10 or more particles in the 
VIRGO and VLS simulations and adopting an equivalent 
mass cut in the GIF simulations (49 or more particles per 
halo in the ACDM and OCDM simulations and 227 or more 
particles in the SCDM and rCDM simulations), such that 
the minimum mass of halos in each simulation is the same 
(this avoids any consequences of possible mass-dependent 
trends in the orbits). 

We find no evidence of any significant difference be- 
tween the velocity distributions constructed from halos with 
10 or more particles and those with 5-20 times more par- 
ticles from the GIF simulations. The measured values of 
X 2 per degree of freedom are scattered around unity and 
are consistent with being drawn from a \ distribution (as 
judged by a K-S test). 

While we would ideally like more extensive tests of the 
effects of particle number 4 we are confident that by using 
halos containing ten or more particles we are obtaining an 
accurate measure of the distributions. 



3.1.4 Radial search limits 

We also wish to test whether our imposed limits on the radial 
separation of halos affects the distributions. To do this we 
use the independent GIF and GIF-11 rCDM z — outputs. 
Velocity distributions are constructed from both simulation 
outputs using radial search limits between Ar = 0.01 and 
Ar = 0.20 in steps of 0.01. We then compute the \ 2 statistic 
comparing the GIF simulation with Ar = 0.20 to the GIF- 
11 simulations with Ar < 0.2, and vice versa. We find that 
the x 2 values stay reasonably constant as the radial search 
limit is decreased, and certainly show no signs of becoming 
significantly larger than unity. As such, we conclude that 
the Ar = 0.2 search limit is sufficiently small to allow an 
accurate determination of the velocity distributions. 



3.2 Trends 

Having established that the techniques employed in this pa- 
per are able to accurately determine the distribution of or- 
bital velocities for infalling satellites we proceed to search 
for any dependence of those distributions on the masses of 
the halos, redshift and cosmology. When testing for such 
dependence we adopt the approach of varying only one vari- 
able at a time, with the hope of isolating the cause of any 
trend we discover. While this is crucial to developing an un- 
derstanding of the trends it significantly limits the number 
of comparisons that we can make. 



4 Ideally we would like a set of simulations identical in all respects 
apart from the number of particles used. This would permit di- 
rect comparisons of the orbital parameters of individual merging 
events to be made. 




V 



Figure 4. Distributions of radial (upper panel) and tangential 
(lower panel) velocities for the GIF and VIRGO SCDM 2 = 
outputs. 

3.2.1 Mass Dependence 

Since our distributions are constructed by combining the 
orbits of all the halos, irrespective of mass, in a given simu- 
lation output it is crucial that we first test for the presence 
of any trends with mass. To test for mass-dependent trends 
we compare the GIF simulations with the VIRGO and VLS 
simulations. These have identical cosmological parameters, 
and we use halos with 10 or more particles in each simula- 
tion. The only difference then is the particle mass and the 
corresponding mass function of dark matter halos. 

We find evidence for mass-dependence in the distribu- 
tions of orbital parameters. Figure |3] shows distributions of 
radial and tangential velocities for GIF and VIRGO SCDM 
models at z = 0. There is a clear difference between the two, 
with the VIRGO simulation showing larger radial and lower 
tangential velocities on average. Unfortunately, our samples 
of mergers remain too small to provide an accurate deter- 
mination of the nature of the mass dependent trends. 

3.2.2 Redshift and cosmology 

We next explore trends with redshift by comparing the 
results of outputs from the same simulation at different 
epochs. Specifically we compute \ 2 f° r pairs of outputs 
which differ by at least 50% in 1 + 2 to ensure that the 
samples are independent. We find strong evidence for differ- 
ences between these samples. However, as the mass function 
of dark matter halos is a function of redshift, we cannot 
disentangle any redshift-dependent trend from the known 
mass-dependent trends. The current simulations do not pos- 
sess enough halos to allow us to select a sub-sample of merg- 
ers by mass at each redshift in order to eliminate this prob- 
lem. We also find significant differences between models with 
different cosmological parameters, but again cannot disen- 
tangle any possible mass-dependent trends. To fully address 
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Figure 5. Distributions of radial (upper panel) and tangential 
(lower panel) velocities for the VIRGO OCDM and SCDM z = 
0.10 outputs. 



Figure 6. Distributions of radial (upper panel) and tangential 
(lower panel) velocities for the VLS and VIRGO ACDM 2 = 
outputs. 



these issues will require a set of custom N-body simulations 
designed to allow us to explore changes in the orbital pa- 
rameter distributions in a controlled manner. (For example, 
the current simulations have a variety of softening lengths, 
which may affect our results. A dedicated set of N-body sim- 
ulations could explore the effects of this parameter on the 
distributions recovered.) 



3.2.3 Group Finding Algorithm 

We test for possible dependence on the group finding algo- 
rithm by comparing distributions of orbital velocities from 
the VLS ACDM simulation with halos found using the FOF 
group finder, to those from the VIRGO ACDM simulation 
with halos found using the SO group finder. We find no 
evidence for any systematic difference between the distribu- 
tions based upon these two group finders, and so use the 
FOF algorithm throughout the remainder of this work. 



3.2.4 Linking Length 

We test for possible dependence on the linking length by 
comparing distributions of orbital velocities from the VLS 
ACDM simulation with halos found using the fixed (varying) 
A, to those from the VIRGO ACDM simulation with ha- 
los found using the varying (fixed) A. The distributions are 
found to be formally inconsistent with one another. Figure|5] 
shows a comparison. With the current statistical precision 
it is difficult to determine the exact nature of the difference 
between fixed and varying A distributions. We will thus not 
explore this further, and will continue to use the varying A 
method. 



3.2.5 Halo Centring Algorithm 

We test for possible dependence on the halo centring algo- 
rithm by comparing distributions of orbital velocities from 
the VLS and VIRGO ACDM simulations with halos found 
using each algorithm (COM and MBP). The distributions 
are again found to be formally inconsistent with one an- 
other. In Figure |S| we show a comparison for the z = 
simulation outputs. The differences between the two distri- 
butions are clearly visible. We find that the COM algorithm 
typically produces distributions of radial and tangential ve- 
locities which peak at lower values than the MBP algorithm. 
As we demonstrated in Fig.0 the COM algorithm can easily 
find an unrealistic origin in both position and velocity space. 
Figure |S] shows that this problem can significantly affect the 
resulting distribution of orbital parameters. We prefer to use 
the more robust MBP algorithm, and do so throughout the 
remainder of this paper. 

3.3 Fitting Functions 

The results presented in this work are potentially of great 
value to any study involving the evolution of the substruc- 
ture population of cold dark matter halos. To facilitate their 
use in this way we provide a simple fitting function which de- 
scribes the two-dimensional distribution of orbital velocities. 
Through simple variable transformations this function also 
describes the distributions of substructure energies, angular 
momenta, eccentricities etc. 

We find that our measured two-dimensional distribu- 
tions of orbital velocities can be reasonably well fit with the 
following fitting function: 

f(vi, v g ) = aive exp [~a 2 {v e - a 9 ) 2 - 6i(ue){u r - b 2 (v e )} 2 ] , (2) 
where 

bi(ye) = a 3 exp [— 04(1*9 - a 5 f] , (3) 
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Figure 7. Distributions of orbital parameters measured in the 
VLS plus VIRGO ACDM z = outputs are shown by crosses. 
Upper left and right-hand panels show distributions of radial and 
tangential velocities respectively. The lower left-hand panel shows 
the distribution of infall angles, while the lower right-hand panel 
shows the two-dimensional distribution of radial and tangential 
velocities (solid contours). Dashed lines show the fitting func- 
tion, while histograms show this function averaged over the same 
bins as used to measure the distributions. The dotted line in the 
lower left-hand panel indicates the distribution of infall angles 
that would occur if correlations between V r and Va were ignored. 



Figure 8. As Figure 171 but for z = 0.5. 



b2{ve) = a 6 exp [-a 7 (ve — a s ) 2 



(4) 
(■>) 



Note that this has a form similar to a two-dimensional 
Maxwell-Boltzmann distribution for the tangential velocity 
and a Gaussian for the radial velocity, as might be expected 
from the results of Vitvitska et al. (2002). However, the 
mean and dispersion of the radial velocity distribution are a 
function of the tangential velocity, as is necessary to account 
for the correlation between these two velocities found in our 
distributions. 

We have fit this function to distributions of orbits taken 
from the combined VLS and VIRGO ACDM simulations 
(the VLS simulation is the only one which provides suffi- 
cient signal to noise to make fitting worthwhile). Figures [7| 
through [5] show distributions of orbital velocities together 
with the fitting function, while Table lists the parameter 
values used in the fits. 



3.4 Other quantities 

Other quantities which characterise the satellite orbits 
are easily derived from the two velocities V T and Ve. For 
convenience, we list below expressions for several other 
orbital parameters in terms of these velocities. 



0.5 1 1.5 2 
V. 





Figure 9. As Figure 171 but for z = 1.0. 



Specific angular momentum: 

J = Ve 

Eccentricity: 



Circularity: 



Specific energy: 



e = V e 



2/2 - Vr - Yl 

2/2-1 



(6) 



(7) 



(8) 



(9) 
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Table 2. Parameters of the fitting function given in eqn. l2l . 
Each column lists parameters which best fit distribution of orbital 
parameters in the combined VLS and VIRGO ACDM simulations 
at the specified redshift. 







Redshift 




Parameter 


0.0 


0.5 


1.0 


a i 


3.90 


4.46 


6.38 


a> 


2.49 


2.98 


2.30 


a.2, 


10.2 


11.0 


18.8 


0,4 


0.684 


1.11 


0.506 


at 


0.354 


0.494 


-0.0934 


ag 


1.08 


1.16 


1.05 


a? 


0.510 


0.261 


0.267 


cix 


0.206 


-0.279 


-0.154 


ag 


0.315 


0.331 


0.157 



Semi-major axis: 

_ h 

a 2/ 2 - V? - V* 
Pericentric distance: 



^"peri — 



Apocentric distance: 



(10) 



(11) 



1 S 



1 - 



(12) 



3.4-1 Eccentricity and semi-major axis 

We have presented results for radial and tangential veloci- 
ties, but of course can just as easily examine invariant pa- 
rameters of the orbits, such as eccentricity and semi-major 
axis. Figure HOI shows distributions of these two parameters 
from the VLS ACDM z = output, together with the dis- 
tributions implied by our fitting function. Our distribution 
of eccentricities is qualitatively, but not quantitatively, in 
agreement with that presented in the first version of the 
preprint (i.e. astro-ph/0309611 version 1, hereafter Khoch- 
far & Burkert (2004) v.l) by Khochfar & Burkert (2004), 
with most orbits being close to parabolic (e = 1). We find 
that almost half of all orbits have e = 1 ± 0.1, a somewhat 
smaller fraction than the 70% given by Khochfar & Burkert 
(2004) v.l. 

In fact, as we show in Fig. 1111 our results are signifi- 
cantly different from those of Khochfar & Burkert (2004) 
v.l. Comparing results from this work with those of Khoch- 
far & Burkert (2004) v.l we find that our results, though 
peaked around e = 1, are more broadly distributed. Khoch- 
far & Burkert (2004) use a different approach to finding 
merging halos than we do and this could potentially in- 

5 Briefly, they locate the progenitors of a given halo at a slightly 
earlier redshift. They then measure the orbital properties of these 
progenitors, providing they are separated by more than the sum 
of their virial radii. To ensure that the apparently merging halos 



fiuence the results obtained. We have implemented Khoch- 
far & Burkert's (2004) methods on the GIF ACDM simu- 
lations to test for any systematic effects caused by the dif- 
ference in methods. We have checked that our implemen- 
tation produces eccentricities identical to theirs (Khochfar, 
private communication). Khochfar & Burkert (2004) v.l 
did not add on the Hubble flow velocity to the motions of 
halos (Khochfar, private communication). Using the Khoch- 
far & Burkert (2004) methods we obtained the distributions 
shown by filled triangles and open squares in Fig. 1111 (Filled 
triangles have no Hubble flow added to halo motions, while 
open squares do have the Hubble flow added.) We find that 
we are able to reproduce the results of Khochfar & Burkert 
when using only halo peculiar velocities in our calculations, 
and are able to reproduce our own results when the Hubble 
flow is included. 

As a second check, we have taken the distribution of or- 
bital circularities found by Tormen (1997), who used tech- 
niques similar to Khochfar & Burkert (2004) , and converted 
these into eccentricities using eqns. © and © (assuming 
fs = 1). Correcting for the fact that orbits with e > 1 are 
not included in the distribution of Tormen (1997) we find an 
eccentricity distribution as shown by the crosses in Fig. 1111 

We conclude that these two different approaches to de- 
termining distributions of halo orbital parameters produce 
consistent results, providing they attempt to measure the 
same quantities. The differences between the distributions 
of eccentricities reported here and by Khochfar & Burkert 
(2004) v.l can be traced to the choice of whether to include 
the Hubble flow in particle velocities (as we did), or to use 
peculiar velocities, as did Khochfar & Burkert (2004) v.l 6 . 



3-4-2 Correlations between pairs of infalls 

We can test for correlations between the infall directions 
of pairs of satellites merging into the same halo. Figure 1121 
shows the distribution of angles <f> between the radius vectors 
of pairs of satellites merging into the same host halo. 7 Note 



are not merely undergoing an unbound "fly-by" they also check 
that the centres of the halos have not moved further apart by a 
later redshift. 

e As a result of discussions regarding these differences, Khochfar 

6 Burkert have revised their calculations to include the Hubble 
flow (see the published Khochfar & Burkert 2004 or version 2 
of the preprint). Their results are then in good agreement with 
those found in this work, as shown by the stars in Fig. 1111 

7 In this and subsequent figures exploring angles between pairs 
of satellites or satellites and the host halo spin we do not in- 
clude our usual weights when constructing the distributions. Our 
weights reflect the fact that, due to the snapshot sampling pro- 
vided by the N-body simulations) we do not see all mergers, but 
only those which occur within a short time after the snapshot. 
When constructing velocity (or eccentricity, semi-major axis etc.) 
distributions, the weighting used corrects for the unobserved pop- 
ulation of mergers. To make the same correction when consider- 
ing the infall angles here we must supplement the weight with 
an assumption about the angular distribution of the unobserved 
mergers. We make the simple assumption that the unobserved 
mergers have the same angular distribution as those which we 
do observe. As such, the resulting angular distribution is found 
from the observed mergers without any weights. Note that this 
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Figure 10. Distributions of eccentricity (left-hand panel) and semi-major axis (right-hand panel) for the VLS plus VIRGO ACDM 2 = 
outputs are shown by the crosses with errorbars. The solid lines indicate the distribution resulting from the fitting formula of eqn. |2] 
The vertical dashed line the left-hand panel indicates parabolic orbits, and so the division between bound (e < 1) and unbound (e > 1) 
orbits. In the right hand panel, negative values of a correspond to unbound orbits. In this case the semi-major axis of the hyperbolic 
orbit is \a\. 



that we have summed the results from all simulation out- 
puts to obtain this distribution. This is permissible as our 
aim here is to search for any deviation from uncorrelated 
infall directions. As such, it does not matter if the differ- 
ent outputs are correlated in different ways — we would still 
see a difference from the null hypothesis of no correlations. 
The distribution appears to differ significantly from that ex- 
pected if there were no correlations between infall directions. 
This correlation between infall directions is qualitatively as 
expected if mergers tend to occur along filaments, i.e. there 
is an enhancement in the number of mergers at small angles, 
£ <J 30° , with a corresponding suppression of mergers with 
angles around 90°. 



3.4-3 Spin alignments 

Finally, we can examine correlations between the infall direc- 
tion and velocity of satellites and the spin angular momen- 
tum vector of the host halo. Figure [T3l shows the resulting 
distributions. We find marginal evidence for deviations from 
a uniform distribution on the sphere. In particular, there is 
a suggestion that merging satellites have a tendency to have 
velocities normal to the spin axis of their host halo. 

To assess the validity of these results will require a bet- 
ter calibration of our errors. For example, the direction of the 
spin vector may be poorly determined for low mass halos, a 
contribution to the errors that we do not take into account. 
(Although this effect should presumably weaken any corre- 
lations, implying that the true correlations are stronger than 
those we measure.) 



assumption may be incorrect — for example, if the angular distri- 
bution correlates with infall velocity — but is at least simple. 



4 DISCUSSION 

We have described methods for determining the orbital pa- 
rameters of dark matter halos at the point of merging with 
a larger system. Previous studies of the orbital properties 
of merging halos have typically considered the orbits after 
merging with the host halo, in which case the orbits will have 
changed due to dynamical friction. Other studies (Tormen 
1997; Khochfar & Burkert 2004) used techniques which are 
restricted to simulations with closely spaced outputs if they 
are to be accurate. Furthermore, we have analysed a sub- 
stantially larger number of orbits than has been previously 
possible to obtain improved statistical precision. This allows 
us to characterise in detail the two-dimensional distribution 
of infall velocities. 

Our analysis pays particular attention to carefully iden- 
tifying halos and their centres. We find that is is important 
to accurately identify the centre of the halo in both posi- 
tion and velocity space, and adopt a similar minimum "en- 
ergy" definition for both of these. We have demonstrated 
that our results are unbiased by effects of particle number or 
radial search limit. In this work, we have focused on the two- 
dimensional distribution of radial and tangential velocities 
which we show has a relatively simple form. A fitting for- 
mula that describes this distribution is presented and should 
prove immensely valuable in future studies of satellite orbits. 

Our methods could be improved upon in several ways. 
A set of simulations run with measurements of orbital pa- 
rameters in mind would allow a better determination of the 
accuracy of our error estimates. More and larger simulations 
would also improve the statistical accuracy of our measure- 
ments and permit us to quantify the trends with, for ex- 
ample, mass that are apparent in the distributions. Finally, 
a more detailed treatment of the evolution of the satellite 
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Figure 13. The distribution of angles between the infall direction (left-hand panel) and infall velocity (right-hand panel) of satellites 
and the angular momentum of the host halo. Points show results measured by summing merger events from all simulation outputs while 
histograms show the expectation when no correlations are present. 
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Figure 11. The distribution of orbital eccentricities. The quan- 
tity shown is the fraction of orbits in each eccentricity bin (i.e. 
following the format of Figure 1 of Khochfar & Burkert 2004). 
Filled circles indicate the results of Khochfar & Burkert (2004) 
v.l, while crosses show the results of Tormen (1997). Open cir- 
cles are results from this work combining all rcdshifts from the 
GIF ACDM simulation using the MBP halo centring algorithm. 
Filled triangles show our implementation of Khochfar & Burkert's 
(2004) methods when no Hubble flow is added to the velocities 
of particles in the N-body simulations, while open squares show 
the same with the Hubble flow added. Stars indicate the results of 
Khochfar & Burkert (2004) which represent the same calculation 
as Khochfar & Burkert (2004) v.l revised to include the Hubble 
flow. 



Figure 12. The distribution of angles between the infall direc- 
tions of pairs of satellites merging onto the same host halo. Points 
show the distribution measured by summing results from all sim- 
ulation outputs while the histogram indicates the expectation for 
uncorrclatcd infall directions. 



orbits (including the effects of an extended, non-spherical 
host halo and dynamical friction) would remove sources of 
systematic error in our measurements. All of these factors 
will be the subject of a future paper. 

We have presented evidence for the presence of trends 
with mass (and, perhaps, with redshift and cosmological pa- 
rameters) in this distribution, although we are currently un- 
able to accurately characterize these trends. Larger samples 
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of N-body halo mergers will allow us to both characterise 
these mass trends and to select sub-samples with a narrow 
range in mass to permit trends with redshift and cosmolog- 
ical parameters to be examined. 

We have also explored the distribution of eccentricities 
and semi-major axes. We find that the eccentricity distribu- 
tion is peaked around parabolic orbits (e = 1). This is qual- 
itatively in agreement with the work of Khochfar & Burk- 
ert (2004) v.l. However, we find quantitative disagreements 
with their distribution of eccentricities. This disagreement 
has been traced to the fact that the Hubble flow was in- 
cluded in our calculations, while it was not included in those 
of Khochfar & Burkert (2004) v.l. Once Hubble flow is in- 
cluded, as in the final version of Khochfar & Burkert (2004) , 
the results of the two studies are in excellent agreement. 
Our distributions of eccentricities and tangential velocities 
are also in good agreement with those from Tormen (1997) 
and Vitvitska et al. (2002) respectively. 

Finally, we searched for correlations between the infall 
directions of pairs of satellites and between the infall posi- 
tions and velocities of satellites and the angular momentum 
of their host halo. We find evidence that satellites infalling 
onto a given host tend to arrive from similar directions, com- 
patible with the hypothesis that (at least some) infall occurs 
along filaments. We find marginal evidence that infall direc- 
tions and direction of motion are aligned with the spin-axis 
of the host halo, although a more thorough study would be 
required to both confirm and interpret this possible correla- 
tion. 

The evolution of sub-structures in cold dark matter ha- 
los is currently a topic of great interest. The tools provided in 
this work should prove of great value in further such studies 
while the techniques described should permit more accurate 
estimates of orbital parameter distributions (including de- 
pendences on halo mass, spin, redshift, cosmology etc.) to 
be constructed. 
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